ELECTROCHEMICAL THIN FILMS AT AND ABOVE THE 
CLASSICAL LIMITING CURRENT 

Q ■ KEVIN T. CHUt AND MARTIN Z. BAZANTt 

fvj 1 Abstract. Wc study a model electrochemical thin film at dc currents exceeding the classical 

' difl'usion-limited value. The mathematical problem involves the steady Poisson-Nernst-Planck equa- 

K ' tions for a binary electrolyte with nonlinear boundary conditions for reaction kinetics and Stern-layer 

^2 capacitance, as well as an integral constraint on the number of anions. At the limiting current, we 

^~i , find a nested boundary layer structure at the cathode, which is required by the reaction boundary 

S jT^ 1 condition. Above the limiting current, a depletion of anions generally characterizes the cathode side 

I ' of the cell. In this regime, we derive leading-order asymptotic approximations for the (i) classical bulk 

' space-charge layer and (ii) another, nested highly charged boundary layer at the cathode. The former 

I I ' involves an exact solution to the Nernst-Planck equations for a single, unscreened ionic species, which 

r~| , may apply more generally to Faradaic conduction through very thin insulating films. By matching 

O ,1 expansions, we derive current-voltage relations well into the space-charge regime. Throughout our 

I ■ analysis, we emphasize the strong influence of the Stern-layer capacitance on cell behavior. 

S: 

fl^ , Key words. Poisson-Nernst-Planck equations, electrochemical systems, limiting current, reac- 

(-H , tion boundary conditions, double-layer capacitance, polarographic curves 
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c/3 . Introduction. Thin-film technologies offer a promising way to construct 

rechargeable micro-batteries, which can be directly integrated into modern electronic 
circuits ^ |21 El El |2] ■ Due to the power-density requirements of many applica- 
tions, such as portable electronics, micro-batteries are likely to be operated under 
at high current density, possibly exceeding diffusion limitation. In a thin film, very 
large electric fields are easily produced by applying only small voltages, due to the 
\^ . small electrode separation, which may be comparable to the Debye screening length. 

[~^ ' Under such conditions, the traditional postulates of macroscopic electrochemical sys- 

^D . tems 13 |S| — bulk electroneutrality and equilibrium double layers — break down 

>0 ' near the classical diffusion- limited current l^l- The mathematical justification for 

~l . these postulates is based on matched asymptotic expansions in the limit of thin dou- 

^\ ' ble layers ^1^1^], which require subtle modifications at large currents. 

"^ . The concept of a "limiting current" , due to the maximum, steady-state flux of 

Q ' diffusion across an electrochemical cell, was introduced by Nernst a century ago |1(J| . 

'^ i but it was eventually realized that the classical theory is flawed, as illustrated in 

^~>' Figure n by numerical solutions to our model problem below. Levich was perhaps 

1-5 ' the first to notice that the assumption of bulk electroneutrality yields approximate 

solutions to the Poisson-Nernst-Planck (PNP) equations, which are not self consistent 
near the limiting current, since the predicted charge density eventually exceeds the 
salt concentration near the cathode (llj . This paradox was first resolved by Smyrl 
and Newman, who showed that the double layer expands at the limiting current, 
C^ ' as the Poisson-Boltzmann approximation of thermal equilibrium breaks down |15| . 

Rubinstein and Shtilman later pointed out that mathematical solutions also exist for 
larger currents, well above the classical limiting value, characterized by a region of 
non-equilibrium "space charge" extending significantly into the neutral bulk |16| . 

The possiblity of supcrlimiting currents has been studied extensively in the differ- 
ent context of bulk liquid electrolytes, where a thin space-charge layer drives nonlinear 






> 



Ph. 



X 



tDcpartmcnt of Mathematics, Massachusetts Institute of Technology, Cambridge, MA 02139 

1 



KEVIN T. CHU AND MARTIN Z. BAZANT 



c ,c 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
X 



Fig. 1. Profiles of the dimensionless potential (top left), electric field (top right), total ionic 
concentration (bottom left), and charge density (bottom right) in three regimes: below the classical 
diffusion-limited current (i = 0.5J, at the limiting current (\ = 1), and above the limiting current 
(i = 1.5j. These are numerical solutions to our model problem with the following dimensionelss 
parameters: e = 0.01, <5 = 0, fcc = 10, i,- = 10. 



electro-osmotic slip. This phenomenon of "electro-osmosis of the second kind" was 
introduced by Dukhin for the nonlinear electrophoresis of ion-selective, conducting col- 
loidal particles JJ_., and Ben and Chang have recently studied it in microfluidics |18| . 
The mathematical analysis of second-kind electro-osmosis using matched asymptotic 
expansions, similar to the approach taken here, was first developed by Rubinstein 
and Zaltzman for related phenomena at electrodialysis membranes |19l EI] . In earlier 
studies, the space-charge layer was also invoked by Bruinsma and Alexander j21| to 
predict hydrodynamic instability during electrodeposition and by Chazalviel "27 in a 
controversial theory of fractal electrochemical growth. 

As in our companion paper on sublimiting currents jH] , here we consider (typically 
solid or gel) thin films, e.g. arising in micro-batteries, which approach the classical 
limiting current without hydrodynamic instability. At micron or smaller length scales, 
the space charge layer need not be "thin" compared to the film thickness, so we also 
analyze currents well above the classical limiting current, apparently for the first 
time. In both regimes, close to and far above the classical limiting current, we derive 
matched asymptotic expansions for the concentration profiles and potential, which 
we compare against numerical solutions. In addition to our focus on superlimiting 
currents and small systems, a notable difference with the literature on second-kind 
electro-osmosis is our use of nonlinear boundary conditions for Faradaic electron- 
transfer reactions, assuming Butler- Volmer kinetics and a compact Stern layer. We 
also analyze the current- voltage relation, thus extending our analogous results for thin 
films below the limiting current [S]. 
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1. Statement of Problem. Before delving into the analysis (and to make this 
paper as self-contained as possible), we review governing equations and boundary 
conditions. We shall focus solely on the dimensionless formulation of the problem. 
For details of the physical assumptions underlying the mathematical model, the reader 
is referred to Ref. ^. 

The transport of cations and anions is described by the steady Nernst-Planck 
equations 

(1) f^ + — (c ^\ =0 

dx'^ dx \ dx J 

while Poisson's equation relates the electric potential to the charge density, 

Here e is a small dimensionless parameter equal to the ratio of the Debye screening 
length to the electrode separation (or film thickness). Note that this formulation 
assumes constant material properties, such as diffusivity, mobility, and dielectric co- 
efficient, and neglects any variations which may occur at large electric fields. The 
factor of 1/2 multiplying the charge density c^ — c_ is present merely for convenience. 
The domain for the system of eqns. Q-® is < a; < 1. 

The two Nernst-Planck equations are easily integrated under the physical con- 
straint that the boundaries are impermeable to anions {i.e. zero flux of anions at 
a; = 0) and taking the nondimensional current density at the electrodes to be 4i: 

/ . X dc+ d(j) 

(4) ^+c+-^=4i 

dx dx 

f5) — - c — = 

dx dx 

Then by introducing the average ion concentration and (half) the charge density 
(6) c=-{c+ + c^) and p=-(c+-c_), 

we can derive a more symmetric form for the coupled Poisson-Nernst-Planck equa- 
tions: 

(7) 

(8) 

(9) 

For this system of one second-order and two first-order differential equations, we 
require four boundary conditions and one integral constraint: 

(10) ^(0)-5e^(0)=0, 

dx 



dc d(f) 
dx dx 


= 2i 


dp dcj) 
dx dx 


= 2i 


^d'^cf) 
' dx^ ■ 


= P- 
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(11) ^(^i) + Se^{l)^v, 

(12) fce[c(0) + p(0)]e"^^('') - ve"""-^'") = i, 

(13) -fcc[c(l) + p(l)]e"=('^(i)-") + ve-""('^(i)-") = 1, 

(14) / [c{x) ~ p{x)]dx = 1, 



The first two boundary conditions, Eqs. H10(I - H11[I . account for the intrinsic capaci- 
tance of the compact part of the electrode-electrolyte interface, as originally envisioned 
by Stern. In these boundary conditions, (5 is a dimensionless parameter which is a 
measure of the strength of the surface capacitance, and v is the total voltage drop 
across the cell. 

The next two boundary conditions, Eqs. ()12|I - H13() . are Butler- Volmer rate equa- 
tions which represent the kinetics of Faradaic electron-transfer reactions at each elec- 
trode, with an Arhhenius dependence on the compact layer voltage. In these equa- 
tions, kc and i^ are dimensionless reaction rate constants and Oc and aa are transfer 
coefficients for the electrode reaction. It is worth noting that ac and aa do not vary 
too much from system to system; typically they have values between and 1 and 
often both take on values near 1/2. 

Finally, the integral constraint, Eq. I|14(l . reflects the fact that the total number 
of anions is fixed, assuming that anions are not allowed to leave the electrolyte by 
Faradaic processes or specific adsorption. It is important to understand that the need 
for an extra boundary condition/constraint refiects that the current-voltage relation- 
ship, i(w), or "polarographic curve", is not given a priori. As usual in one-dimensional 
problems 0, it is easier to assume galvanostatic forcing at fixed curent, i, and then 
solve for the cell voltage, w(i), by applying the boundary condition ljll() . rather than 
the more common case of potcntiostatic forcing at fixed voltage, v. For this reason, 
we take the former approach in our analysis. For steady-state problems, the two kinds 
of forcing are equivalent and yield the same (invertible) polaragraphic curve, i(u) or 

For some of our analysis, it will be convenient to further simplify the problem by 
introducing the dimensionless electric field, E = —-^. This transformation is useful 
because three of the five independent constraints can be expressed in terms of these 
variables, without explicit dependence on 0(a;), namely the two Butler- Volmer rate 
equations, 

(15) fc,(c(0) + p(0))e^"=*^^(") - ve""*^^(°) = i, 

(16) -fcc(c(l) -I- p(l))e"=^^^(i) + i,e-""*^^(i) = 1, 

and the integral constaint on the total number of anions, Eq. (|14|) . The potential is 
recovered by integrating the electric field and applying the Stern boundary conditions 
Eqs. (P|l and ifTT)) . 

2. Unified Analysis at All Currents. 

2.1. Master Equation for the Electrostatic Potential. We begin our anal- 
ysis by reducing the governing equations, Eqs. 10 through (|5J), to a single master 
equation for the electrostatic potential. Substituting Eq. ^ into Eq. iQ and inte- 
grating, we obtain an expression for the average concentration 

(17) c(x)=c„ + 2ix+^(^g 
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Then by applying the integral constraint, Eq. ()14|l. we find that the integration con- 
stant, Cg , is given by 



(18) 



Co = (l-i)-e' 



dx 



dx 



dx 



dx 



Note that when the electric field is 0(1), Eqns. H17() and 1)18(1 reduce to the leading- 
order concentration in the bulk when i is sufficiently below the limiting current ^. 
We can now eliminate p and c from Eq. (jH]) to arrive at a single master equation for (p 



+ (e, + 2..)^.2. 



(19) 


d^ct) 1 fd(f>\ 
dx^ 2 \dx J 


3" 


or equivalently for the electric field E 


(20) 


e' 


\d'E 11 
dx^ 2 


- 



(Cq -I- 2ia;) E = 2i. 



Once this equation is solved, the concentration, c, and charge density, p, are computed 
using Eq. ((T7|l and Poisson's equation, Eq. ©. 

The master equation has been derived in various equilivalent forms since the 
1960s. Grafov and Chernenko |53] first combined Eqs. Q, (0) and Q to obtain a 
single nonlinear differential equation for the anion concentration, c_, whose general 
solution they expressed in terms of Painleve's transcendents. The master equation for 
the electric field, Eq. (|20|l . was first derived Smyrl and Newman ^^, in the special 
case of the classical limiting current, where i = 1 and c^ = 0, where they discovered 
a non-equilibrium double layer of width, e^'^, which is apparent from the form of 
the master equation. We shall study the general electric-field and potential equations 
for an arbitrary current, i, focusing on boundary-layer structure in the limiting and 
super limiting regimes. 

2.2. EfRcient Numerical Solution. To solve the master equation for the elec- 
tric field with the boundary conditions and integral constraint, we use the Newton- 
Kantorovich method 25 . Specifically, we use a Chebyshev pseudospectral discretiza- 
tion to solve the linearized boundary- value problem at each iteration |25[ 1^ . Our 
decision to use this method is motivated by its natural ability to resolve boundary 
layers and its efficient use of grid points. We are able to get accurate results for many 
parameter regimes very quickly (typically less than a few minutes on a workstation) 
with only a few hundred grid points, which would not be possible at large currents 
and/or thin double layers using a naive finite-difference scheme. It is important to 
stress that the boundary conditions and the integral constraint are explicitly included 
as part of the Newton-Kantorovich iteration. Therefore, the linear BVP solved in 
each iteration is actually an integro- differential differential equation with boundary 
conditions that are integro-algebraic equations. 

To ensure convergence at high currents, we use continuation in the current density 
parameter, i, and start with a sufficiently low initial i that the bulk electroneutral 
solution is an acceptable initial guess; often, initial i values relatively high compared 
to the diffusion-limited current are acceptable. Continuation in the 5 parameter is 
also sometimes necessary to compute solutions at high 5 values. 

The results of the numerical method are presented in the figures below and in 
Ref. [5] to test our analytical approximations obtained by asymptotic analysis. 
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2.3. Recovery of Classical Results Below the Limiting Current, i <C 

1 — e^' ■^. In the low-current regime, the master equation admits the two distinguished 
Umits around x = that arise in the classical analysis: x = 0(1) and x = 0(e). 
When X = 0(1), we find the usual bulk electric field from Eq. (|19|l and the bulk 
concentration from Eq. I|17() . When x — 0(e), the master equation can be rescaled 
using x — ey to obtain 

d^cj) 1 /d(j)Y d(j) ^ d(j) ^ 

which is equivalent to the classical theory at leading order [Hj. In particular, the 
Gouy-Chapman structure of the double-layer can be derived directly from the Smyrl- 
Newman equation in this limit |23| . 

The anode boundary layer comes from a similar 0(e) scaling around x — 1. Note 
that in the i <C 1 — e^/^ regime, the scaling x = 0(e^/'^) is not a distinguished limit 
because the c^ ( ^ j term would dominate all other terms in Eq. (|19|l . 

3. Nested Boundary Layers at the Limiting Current, i = 1 0(e^/'^). In 

this section, we show that a nontrivial nested boundary-layer structure emerges at 
the classical limiting current when general boundary conditions are considered. 

3.1. Expansion of the Double Layer Out of Equilibrium. As discussed in 
the companion paper [Hj , the classical analysis breaks down as the current approaches 
the diffusion- limited current. Mathematically, the breakdown occurs because a new 
distinguished limit for the master equation appears as i — > 1. Rescaling the master 
equation using x = e^'^z gives us 

^ ' dz^^ 2\dz) + £2/3 dz + dz " ' 

which implies that we have a meaningful distinguished limit if c^ = 0{€^'^) or, equiv- 
alently, i = 1 — 0{e^/^). As first noticed by Smyrl and Newman ^2|j the appearance 
of the t^'^ scaling corresponds to the expansion of the diffuse charge layer that arises 
at the classical limiting current (see Figure |21). In this regime, the double layer is no 
longer in Poisson-Boltzmann equilibrium at leading order. 

Unfortunately, at this scale, a/Uerms in Eq. H22|l are 0(1), so we are forced to solve 
the full equation. Although general solutions can be expressed in terms of Painleve's 
transcendents '8, 18, 24 , these are not convenient for applying our nonlinear boundary 
conditions or obtaining physical insight. Even when c^ = o{e^/'^), we are left with a 
complicated differential equation which does not admit a simple analytical solution. 
However, in the case c^ — o{e^'^), it is possible to study the asymptotic behavior 
of the solution in the limits z — > and z — > oo by considering the behavior of the 
neighboring asymptotic layers. 

3.2. Nested Boundary Layers when |1 — 1| = o{e'^^'^). The appearance of the 
new distinguished limit for i = 1 — 0(6^/^) does not destroy the ones that exist in the 
classical analysis. In particular, the 0(e) boundary layer at a: = does not vanish. 
This inner layer was overlooked by Smyrl and Newman because they assumed a fixed 
surface charge density given by the equilibrium zeta potential 15 , rather than more 
realistic boundary conditions allowing for surface-charge variations. 

In the general case, a set of nested boundary layers exists when the current is near 
the classical limiting current. For convenience, we shall refer to the x — 0{e'^''^) and 
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Fig. 2. Numerical solutions for the dimensionless electric field E{x) at current densities 
of I = 0.9 and i = 1.0 demonstrating the expansion of the diffuse layer at the limiting current 
(kc = 1, ir = 2, (5 = 0.1 and e = 0.0001/ For reference, the vertical line shows where x - 



,2/3 



the X — 0{e) regions as the "Smyii-Newman" and "inner diffuse" layers, respectively. 
It is important to realize that without the 0(e) layer, it would be impossible, in 
general, to satisfy the boundary conditions of the original problem. To see this, 
consider the reaction boundary condition at x = 0, Eq. (|12|l . To estimate the c and p 
at the electrode surface, we rescale Eq. H17|l and Poisson's equation using x 
to obtain 



g2/3^ 



(23) 

(24) 



2ie2/3 



,2/3 



2 + 



dz 



^2/3 : 



dz2 



which means that the concentration and charge density are both Oi^e^/^) since c^ = 
o{e^/^) when |1 — i| = o{e^/^). Then, from the Stern boundary condition, we have 
0(0) = —deE = —5e^/'^E = 0((5e^/'^). Plugging these estimates into the reaction 
boundary condition, we find 



(25) 



fc,O(e2/3)ga.5.V3^(0) 



IrC 



^a„<5e^/^£;(0) 



o(i). 



For this equation to be satisfied, 5 is constrained to be huge for e near 0: 5 = 
O (jloge^/'^l /e^^^) ^ 1. But 5 is an independently assigned physical parameter, so 
there is no reason that it should have to satisfy this constraint. Thus, we are lead to 
the conclusion that there must be a nested inner boundary layer in order to satisfy 
the Butler- Volmer boundary condition when 5 is small compared to I log e^/^ I /e^/^ . In 
other words, for insufficiently large 5, the reaction boundary condition requires that 
the cation concentration is 0(1) at x = 0, which implies the existence of a boundary 
layer between the Smyrl-Newman layer and the boundary. As can be seen in figure 
|31 the cathode layer concentration decreases significantly as 5 is increased. 

To analyze Eq. H22|) . it is convenient to focus on the electric field rather than 



the potential, 
becomes 

(26) 



In terms of the scaled electric field, E{z) 



d(l> 



e^/^E{x), Eq. ^ 



(Pe 1 



fi 



dz^ 



- -E^ -2i zE 



l'-#^ 
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which we shall refer to as the "Smyrl-Newman equation". From Eq. (71) in Ref. [^, 
we know that the first few terms in the expansion for the bulk electric field at the 
limiting current are 

- 1 3e2 llle4 6045e6 

1/1 3 111 6045 



^^''^ £2/3 Vz ' 4z4 ' I6z7 ' 32z 

Since the second series is asymptotic for z ^ 1, the expansion in the bulk is valid 
for X ^ e^'^. In order to match the solution in the Smyrl-Newman layer to the 
bulk, we expect the asymptotic solution to Eq. H25|l as 2; ^ oo to be given by the 
expression in parentheses in Eq. lf?7|l . We could also have arrived at this result by 
directly substituting an asymptotic expansion in 1/z and matching coefficients. As 
we can see in Figure the leading order term in Eq. 127|l is a good approximation to 
the exact solution in the bulk and is matched by the solution in the Smyrl-Newman 
layer as it extends into the bulk. 

We now turn our attention towards the "inner diffuse" layer which gives us the 
asymptotic behavior of the Smyrl-Newman equation in the limit z — > 0. Introducing 
the scaled variables y — x/e = z/e^f^ and E = eE = e^/^E, Eq. (|26|l becomes 

(28) 0-^^'-2.(yi^ + l).C.i?. 

Near the limiting current {i.e. c^ — 0{e'^^'^)), E satisfies ^-f = ^& at leading order 

with the boundary condition i^J — > as y — > 00 from the matching condition that E 
remains bounded as z — > 0. Integrating this equation twice with the observation that 
^ > gives us 

2 



(29) E{y) ^ - 

y -r u 

where & is a constant determined by applying the Butler- Volmer reaction boundary 
condition at the cathode. We can estimate c{y) and p{y) by substituing Eq. H29|l into 
Eq. lfT7|) and Poisson's equation to find 



LEixf^c, + 2iey+^E'-'' 



(30) c(y) ^c, + 2ix + -E{xY = c„ + 2xey + -E{yY = -3 + 0(e) 

^ [y + o) 

^, , ^dE dE 2 

31 p{y)^e''—^ — ^- -2+Oe. 

dx dy [y + h) 



Therefore, h satisfies the following transcendental equation at leading order: 

(32) K^e^""=^'^ = 1 + ve-2"<.<5/b. 

While this equation does not admit a simple closed form solution, we can compute 
approximate solutions in the limits of small and large 5 values. In the small 5 limit, 
we can linearize Ea. (|32|l and expand fo in a power series in 5 to obtain 



(33) b^2J-^^ + 5[a, + ^^^\+0{S'). 
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Fig. 3. Numerical solutions (solid lines) for the dimensionless electric field E(x) and concen- 
tration c{x) at the classical diffusion- limited current (i = 1) compared with leading order asymp- 
totic approximations (dashed and dot-dashed lines) for fcc = 1, ir = 2, e = 0.01 and S = 0.1, 10. 
The leading order bulk approximations for E(x) and c{x) are given by Eq. {2Tf and c(x) = 2ix, 
respectively. In the diffuse layer, the leading order approximations are given by Eqs. 1296 and 
lyut . For the S = 10 curves, the difference between the dashed and dot-dashed curves is that the 
dahsed curve uses an approximate value for B given by Eq. I35{l while the dot-dashed curve uses 
a B value calculated by numerically solving Eq. ISiX . For reference, the vertical lines show where 
X = e and x = e^ " . The thin anode diffuse layer field is not shown. 



At the other extreme, for (5^1, Eq. 1321) can be approximated by 



(34) 



^ ^ 2a.S/b _ 



Then, using fixed-point iteration on the approximate equation, we find that 

2ar5 



(35) 



log K — 2 log log K + O (log log log (5^) 



where k = la^iJ^/fcc- Figure |21 shows that the leading order approximation Eq. (|29|l is 
very good in the inner diffuse layer as long as an accurate estimate for h is used. While 
the small 5 approximation for b is amazingly good (the asymptotic and numerical 
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solutions are nearly indistinguishable), the large S estimate for b is not as good but 
is only off by an 0(1) multiplicative factor. 

Before moving on, it is worth noting that the asymptotic behavior of the con- 
centration and charge density in the Smyrl-Newman layer as z — > and z —f oo 
suggest that the charge density is low throughout the entire Smyrl-Newman layer. 
Figure 131 shows how the Smyrl-Newman layer acts as a transition layer allowing the 
bulk concentration to become small near the cathode while still ensuring a sufficiently 
high cation concentration at the cathode surface to satisfy the reaction boundary 
conditions. The transition nature of the Smyrl-Newman layer becomes even more 
pronounced for smaller values of e. 

4. Bulk Space Charge Above the Limiting Current, 1 + 0(e^/'^) <C i ^ 
0(l/e). As current exceeds the classical limiting value, the overlap region between 
the inner diffuse and Smyrl-Newman layers grows to become a layer having 0(1) 
width. Following other authors |16[I22| . we shall refer to this new layer as the ''space- 
charge" layer because, as we shall see, it has a non-negligible charge density compared 
to the rest of the bulk. Therefore, in this current regime, the central region of the 
electrochemical cell is split into two pieces having 0(1) width separated by a o(l) 
transition layer. 

In the bulk, the solution remains unchanged except that c^ cannot be approxi- 
mated by 1 — i; the contribution from the integral term is no longer negligible. The 
need for this correction arises from the high electric fields required to drive current 
through the electrically charged space-charge layer. With this minor modification, we 
find that the bulk solution is 

c{x) = c^ + 2ia; 
(36) E{x) = — ^ — 

where Xo = — Co/2i is the point where the bulk concentration vanishes (see Figure ^J. 
In between the two 0(1) layers, there is a small transition layer. Rescaling the 
master equation using the change of variables z = (x — Xo)l(^l'^ and E{z) = e^l'^E(x\ 
we again obtain the Smyrl-Newman equation, Eq. (|26|l . with right hand side equal 
to zero. As before, we find that the solution in the transition layer approaches — l/z 
as z — > oo. In the other direction as z — > — cx), we will find that the appropriate 
boundary condition \b E ^> —2\J\\z\ to match the electric field in space-charge layer. 

4.1. Structure of the Space-Charge Layer. Physically, we could argue that 
the concentration of ions in the space-charge layer is very small {i.e. zero at leading 
order) because the layer is essentially the result of stretching the ionic content of the 
overlap between the inner diffuse and Smyrl-Newman layers, which is small to begin 
with, over an 0(1) region. This physical intuition is confirmed by the numerical 
solutions shown in Figures 01 El and Therefore, using Eq. H17|l . we obtain the 
leading order solution for the electric field 

— 2\/l iXn — x\ 

(37) E y_^_^ ^_ 

e 
Note that the magnitude of the field is exactly what is required to make the integral 
term in c^ an 0(1) contribution. From this formula, it is easy to compute the charge 
density in the space-charge layer 



^dE 
(^«) '' = '^x 
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Fig. 4. Numerical solutions (solid lines) for the dimensionless electric field E{x), average 
concentration c{x), and charge density p(x) above the diffusion-limited current (\ = 1.5j compared 
with leading order asymptotic approximations (dashed lines) for fcc = 1, ir = 2, e = 0.01, and 
& = 0.1,10. The leading order bulk approximations are given by Eqs. lS6Xl . In the space-charge 
layer, the leading order electric field is given by Eq. X'Jl^ , and leading order concentration is 0. 
Finally, Eqs. |57| ) and 15 fX are the diffuse layer asymptotic approximations for the electric field 
and concentration, respectively. For reference, the vertical lines show where x = e and x = Xo. 



which is an order of magnitude larger than the O(e^) charge density in the bulk. The 
0(e) charge density also implies that the concentration must be at least 0(e) because 
the anion concentration, c — p, is positive. 

With the electric field given by Eq. (|37|l , we can determine the values of Xo and c^ 
by solving the system of equations given by the definition of Xq and c^. Using Eq. H18|l 
to calculate c^ and noticing that the leading order contribution to the integral comes 
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Fig. 5. Numerical solutions for the dimensionless cation and anion concentrations above 
the diffusion-limited current (\ = l.b) for fcc = 1, ir = 2, e = 0.01, and 5 = 0.1, 10. For reference, 
the vertical lines show where x = e and x = Xo- 



from the space-charge layer, we obtain 

(39) c,^l-i{l + xl). 

Combining this result with Xo = — Cq/2i, we find that 



(40) 



1-1-1/2 



(.--.), 



which can be substituted into Eqs. (|36|l and (|37|l to yield the leading order solutions 
in the bulk and space-charge layers. It should be noted that the expression for Xo is 
consistent with the estimate for the width found by Bruinsma and Alexander |23 and 
Chazaviel |35| in the limits i— 1 ^ 1 and small space-charge layer {xo <C 1), although 
our analysis also applies to much larger voltages. 

The results obtained via physical arguments in the previous few paragraphs moti- 
vate an asymptotic series expansion for E whose leading order term is 0(l/e). More- 
over, because we want to be able to balance the current density at second-order, we 
choose the second-order term to be 0(i). Thus, we have 



(41) 



E 



1 



--E/-1 + Eqi - 
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Fig. 6. Numerical solutions (solid lines) for the dimensionless electric field E{x), average 
concentration c{x), and charge density p{x) far above the diffusion-limited current (\ = 10. Oj 
compared with leading order asymptotic approximations (dashed lines) for kc = 1, ij- = 2, e = 0.01, 
and 5 = 0.1. Each field is shown twice: (1) with x on log scale to focus on the cathode region 
and (2) with x on a linear scale to emphasize the interior of the cell. Note that it = 0.1, so 
the asymptotic approximations are not as good as at lower current densities. For reference, the 
vertical lines show where x = e and x = Xo. 



Note that in this asymptotic series, the first term only dominates the second term as 
long as 1 ■€. 1/e, so the following analysis only holds for current densities far below 
0(l/e). Figure El illustrates the breakdown of the leading-order asymptotic solutions 
at very high current densities. While the qualitative features of the asymptotic ap- 
proximation are correct {e.g. shape of E{x) in the diffuse layer and slope of c{x) in 
the bulk), the quality of the approximation is clearly less than at lower values of i. 

The key advantage of a more systematic asymptotic analysis is that we are able to 
calculate the leading-order behavior of the space-charge layer concentration c, which 
is not possible with only knowledge of the leading order behavior for the electric field. 
Substituting Eq. H41|) into the master equation H2U|) , it is straightforward to obtain 



(42) 



2 /— 1 
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Using this expression in Eq. I|17[l. we find the dominant contribution to c is exactly 
the same as p: 



(43) 



Since c_ = c — p, this result leads to an important physical conclusion — The space- 
charge layer is essentially depleted of anions, c_ = o(e), as is clearly seen in Figures^] 
and |S1 This contradicts our macroscopic intuition about electrolytes, but, in very 
thin films, complete anion deplection might occur. For example, in a microbattery 
developed for on-chip power sources using the Li/Si02/Si system, lithium ion conduc- 
tion has recently been demonstrated in nano-scale films of silicon oxide, where there 
should not be any counter ions or excess electrons (HI- 

At leading order as e ^ 0, the anion concentration, c_, can be set to zero in the 
space-charge layer, leaving the following two governing equations: 



(44) 
(45) 



dx dx 

2^20 1 



As with binary electrolyte case, these equations can be reduced to a single equation 
for the electric potential: 



(46) 



d^(f) d'^ct) d(f) 

dx^ dx^ dx 



2i 
72 ■ 



Integrating this equation once, we obtain a Riccati equation for -^ 

d'^<t> 1 fd4)V 2i, , , 

(4^) d^+2UJ =-^(x -.„) + /., 

where h is an integration constant. Using the transformations 



(48) 



= 0/2 



we find that u satifies Airy's equation 

'd^ 



(49) 



'^2/3 '^ ^°' + 2i2/3 ' 



ZU = 0. 



Thus, the general solution for (/)(a;) is 
(50) </)(a;) = 2 log 



jl/3 \ /jl/3 

aiAi ( -— {xo - x) + I3h\ + a2Bi ( -— {xo - x) 



Ph 



where ai and 02 are constants determined by boundary conditions and P = ^2/3 ■ 

To simplify this expression, note that in the limit e ^ 0, the potential drop 
between x = Xq and a; = is approximately 



(51) 



(xo)- 0(0) -2 log 



aiAi{0) + a2Bi{0) 



aiAi (^^) + a2Bi (^^) 
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Now, using the large argument behavior of the Airy functions, we see that as e ^ 0, 
the argument of the logarithm approaches zero. Thus, we are lead to the conclusion 
that the electric potential at x = Xq is less than at a; = 0. But, this is completely 
inconsistent with our physical intuition and the numerical results, which show that 
4>(xo) — 0(0) > 0. Therefore, it must be the case that 02 ~ so that 

(52) 0(a;) = 2 log aiAi l^j^{xo-x)+ f3h 

and 

2ii/3 Ai' (^^ (xo -x)+ I3h 



(53) E{x) 



.2/3 



A-i [^ {xo - x)+ f3h 



Finally, by using the asymptotic form of Ai(z) and Ai'{z) as z — > cx), we find that in 
the e — > limit, the leading order approximation for the electric field when the region 
is depleted of anions is exactly Eq. H37|l . 

That equivalence of the single-ion equations and the full governing equations 
at leading-order mathematically confirms the physically interpretation of the space- 
charge layer as a region of anion depletion. From an alternative perspective, it also 
reminds us that under extreme conditions, it may be necessary to rethink our assump- 
tions about what physical effects are dominant. 

4.2. Boundary Layers Above the Limiting Current. To complete our anal- 
ysis of the high-current regime, 1 + ©(e^/"^) <C 1 ^ 0(l/e), we must consider the 
boundary layers. At the anode, all fields are 0(1), so we recover the usual Gouy- 
Chapman solution with the minor modification that ci = 2^/1 which is the value c 
takes as X —f 1. The cathode structure, however, is much more interesting because 
it is depleted of anions (see Figure [SJ. To our knowledge, this non-equilibrium inner 
boundary layer on the space-charge region, related to the reaction boundary condition 
at the cathode, has not been analyzed before. 

As in the space-charge layer, the leading-order governing equations in this layer 
are those of a single ionic species with no counterions Eqns. (|44|1 and H45|) . Rescaling 
those equations using x — ey, we obtain 

dc I d<i) 

(54) ^+c+--=4iew0 

ay ay 



(55) -3^=0^+- 



(f4> _ 1 

From these equations, it is immediately clear that the cations have a Boltzmann 
equilibrium profile at leading order: c+ cx e~'^^'^'. As in the analysis for the space- 
charge layer, it is possible to find a general solution to Eqns. H54|l and H55(l . By 
combining these equations and integrating, we find that the potential in the cathode 
boundary layer has the form 

(56) (j) ~ log [sinh^ {py + q)] + r, 

where p, q, and r are integration constants. Therefore, the electric field and concen- 
tration are 

(57) E{y)^~2pcoth{py + q) 

(58) c{y) = lc+{y) ^ —^ 



2 sinh (jiy -|- q) 
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Matching the electric fields in the diffuse and space-charge layers, we find that 
p ^ y/^xZ- Note that because p = 0(-\/T), the electric field in the diffuse charge layer 
is 0{^/i/e) which is same order of magnitude as in the space-charge layer. To solve 
for q, we use the expression for p in the cathode Stern and Butler- Volmer boundary 
conditions, which leads to the following nonlinear equation: 

(59) -; — ^— e:xJp(2acS^/TxZcothq) — irC:iq){—2aaS^/TxZcothq) — i. 



sinh q 

In the limit of small S, we can use fixed-point iteration to obtain an approximate 
solution 



, , 1 / / kciXo exp (2acS^/iXo cothoo) 

(60) q - sinh"^ 1^ I c o yy c ^/ o Ho) 



1 -I- ir exp (— 2aQ(5y^i2?^cothqo) / 

where qo has the same form as q with (cothoo) set equal to 1. For (5^1, the leading 
order equation is 

(61) 2~ 6xp {2ac5^\Xo cothq) ~ i, 

sinh q 

which implies that g ^ 1 so that the left-hand side can be small enough to balance 
the current. Thus, by using cothq « 1 and sinhg « exp{q)/2, we find that q ~ 
ac5^J\Xo + 5 log(16fcca;o)- The agreement of these asymptotic approximations with 
the exact solutions in the diffuse charge layer is illustrated in Figure 0] 

5. Polarographic Curves. We are now in a position to compute the leading- 
order behavior of the polarographic curve at and above the classical limiting current. 
Recall that the formula for the cell voltage is given by 

(62) V = -SeEiO) + f -E{x)dx - 5eE{l). 

Jo 

The integral is the voltage drop through the interior of the cell and the first and last 
terms account for the potential drop across the Stern layers. 

At the limiting current, i = 1, we can estimate the voltage drop across the cell by 
using the bulk and diffuse layer electric field to approximate the field in the Smyrl- 
Newman transition layer to obtain 

(63) V ~ -SeE{0) + [ -E{x)dx + [ -~E{x)dx - 5eE{l) 

Jo Je2/3 

(64) ^2--F21og(^ j-gloge. 

Notice that in the small S limit, this expression reduces to u ~ — g Ine as e — > which 
can be observed in the polarographic curves in Figure 6 in [Hj. Table Q compares this 
approximation with the exact cell voltage for a few e and S values. For small e values 
(e < 0.01), the asymptotic approximations are fairly good (within 5% to 10%). 

Above the limiting current, the space-charge layer makes the dominant contribu- 
tion to the cell voltage. Using Eqns. (|36|l and (|37|l in the formula for the cell voltage, 
we find that 

(65) V ~ ^ (l - 1-1/2^^/' + 2,5 (i - ^/I)'^" cothg - i logi - 2/31oge. 
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Table 1 
Comparison of the asymptotic approximations Eqs. |6'4| ) and I6'5I ) with numerically calculated 
values for the cell voltage at various e and 5 values. These cell voltages were computed with fee = 1 
and ir = 2. 





1 


= l.U 




e 


5 


'^exact 


'^asym 


le-4 


0.01 


13.125 


12.101 


le-4 


1.00 


13.222 


12.374 


le-4 


10.0 


14.290 


13.571 


le-3 


0.01 


10.165 


9.146 


le-3 


1.00 


10.277 


9.475 


le-3 


10.0 


11.552 


10.890 


le-2 


0.01 


7.339 


6.303 


le-2 


1.00 


7.479 


6.729 


le-2 


10.0 


9.228 


8.465 


le-1 


0.01 


4.922 


3.649 


le-1 


1.00 


5.005 


4.219 


le-1 


10.0 


7.995 


6.327 



1 = 1.5 



e 


5 


'^exact 


'^asym 


lc-4 


0.01 


1297.799 


1289.621 


le-4 


1.00 


1297.048 


1291.101 


le-4 


10.0 


1305.318 


1300.129 


le-3 


0.01 


140.207 


132.790 


le-3 


1.00 


139.450 


134.270 


le-3 


10.0 


147.717 


143.299 


le-2 


0.01 


22.434 


15.725 


le-2 


1.00 


21.624 


17.206 


le-2 


10.0 


29.886 


26.234 


le-1 


0.01 


9.479 


2.637 


le-1 


1.00 


7.790 


4.118 


le-1 


10.0 


16.088 


13.146 



The first two terms in tliis expression estimate the vohage drop across the space-charge 
and the cathode Stern layers, respectively. The last two terms are the subdominant 
contribution from the bulk where we have somewhat arbitrarily taken x = Xo + e^'^ as 
the boundary between the bulk layer and the Smyrl-Newman transition layer. Notice 
that we ignore the contribution from the cathode diffuse and Smyrl-Newman layers. 
It is safe to neglect the diffuse layer because it is an 0(1) contribution. However, the 
Smyrl-Newman layer has a non-neglible potential drop that we have to accept as error 
since we do not have an analytic form for the solution in that region. 

Figure \7\ shows that the asymptotic polarographic curves are quite accurate for 
sufficiently small e values. In Tabled we compare the results predicted by the asymp- 
totic formula with numerical results for a few specific values of e and 6. It is interesting 
that the approximation is also better for large S values (we explain this observation 
in the next section). Also, while the loge term is subdominant, it makes a significant 
contribution to the cell voltage for e values as small as 0.01. 

As with the width of the space-charge layer, Xo, our expression for the cell volt- 
age, Eq. (|65|l . is consistent with the results of Bruinsma and Alexander [21] and 
Chazaviel 22 , near the limiting curent, i -^ 1+, while remaining valid at much larger 
currents, i = 0(l/e) 

6. Effects of the Stern-Layer Capacitance. The inclusion of the Stern layer 
in the boundary conditions allows us to explore the effects of the intrinsic surface 
capacitance on the structure of the cell. From Figures 13 through [SJ we can see that 
smaller Stern-layer capacitances (i.e. larger d values) decrease the concentration and 
electric field strength in the cathode diffuse layer. This behavior arises primarily from 
the influence of the electric field on the chemical kinetics at the electrode surfaces. 
When the capacitance of the Stern layer is low, small electric fields at the cathode 
surface translate into large potential drops across the Stern layer, Eq. (|10|l . which 
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Fig. 7. Comparison of numerical polarographic curves (dashed lines) with leading-order 
asymptotic approximations (solid lines) given in Eq. I6'5I ) for several values of e with <5 = 1.0, 
fcc = 1 and \r = 2. For e = 0.001, the numerical and asymptotic polarographic curves are indis- 
tinguishable on this graph. For reference, the vertical, dashed line shows the classical diffusion- 
limited current i = 1 . 



help drive the deposition reaction, Eq. H12() . As a result, neither the electric field nor 
the cation concentration need to be very large at the cathode to support high current 
densities. These results confirm our physical intuition that it is only important to 
pay attention to the diffuse layer when the Stern layer potential drop is negligible (i.e 
6 « 1). 

At high currents, another important effect of the Stern layer capacitance is that 
the total cell voltage becomes dominated by the potential drop across the Stern layer 
at large S values {i.e. small capacitances). This behavior is clearly illustrated in 
Figure IHl Notice for currents below the classical diffusion- limited current, the total 
cell voltage does not show a strong dependence on S. However, for i > 1, the total 
cell voltage increases with 6 - the increase being driven by the strong S dependence 
of the Stern voltage. 

7. Conclusion. In summary, we have studied classical problem of direct current 
in an electrochemical cell, focusing on the exotic regime of high current densities. A 
notable new feature of our study is the use of nonlinear Butler- Volmer and Stern 
boundary conditions to model a thin film passing a Faradaic current, as in a micro- 
battery. We have derived leading-order approximations for the fields at and above 
the classical, diffusion-limited current, paying special attention to the structure of 
the cathodic boundary layer, which must be present to satisfy the reaction boundary 
conditions. In our analysis of superlimiting current, we have shown that the key 
feature of the bulk space-charge layer is the depletion of anions. Our exact solution of 
the leading-order problem in the space-charge region, Eq. H50() . could thus also have 
relevance for Faradaic conduction through very thin insulating films. 

Using the asymptotic approximations to the fields, we are able to derive a current- 
voltage relation, Eq. H65|l . which compares well with numerical results, far beyond the 
limiting current. Combined with the analogous formulae in the companion paper [^, 
which hold below the limiting current, we have essentially analyzed the full range of 
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Fig. 8. These graphs break the total cell voltage into contributions from the cell interior and 
the Stem layer as a function of 5 for e = 0.01, kc = 2, and ir = 2. Note that at and above the 
classical limiting current, the Stern layer voltage dominates the total cell voltage for large values 
of 5. 



the current-voltage relation. These results could be useful in interpreting experimental 
data, e.g. on the internal resistance of thin-film niicrobatteries. 

A general conclusion of this study is that boundary conditions strongly affect 
the solution. For example, the Stern-layer capacitance, often ignored in theoretical 
analysis, plays an important role in determining the qualitative structure of the cell 
near the cathode, as well as the total cell voltage. The nonlinear boundary condi- 
tions for Butler- Volmer reaction kinetics also profoundly affect charge distribution 
and current-voltage relation, compared to the ubiquitous case of Dirichlet bound- 
ary conditions. The latter rely on the asumption of surface equilibrium, which is of 
questionable validity at very large currents. 

We leave the reader with a word of caution. The results presented here are valid 
mathematical solutions of standard model equations, but their physical relevance 
should be met with some skepticism under extreme conditions, such as superlimiting 
current. For example, the PNP equations are meant to describe infinitely dilute 
solutions in relatively small electric fields 1271 I28| . Even for quasi-equilibrium 
double layers, their validity is not so clear when the zeta potential greatly exceeds the 
thermal voltage, because co-ion concentrations may exceed the physical limit required 
by discreteness (accounting also for solvation shells) and counter-ion concentations 
may become small enough to violate the continuum assumption. Large electric fields 
can cause the permittivity to vary, by some estimates up to a factor ten, as solvent 
dipoles become aligned. Including such effects, however, introduces further ad hoc 
parameters into the model, which may be difficult to infer from experimental data. 
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